######  Observable Bounds of Rationality and Credibility, Journal of Politics #######

###### Analysis Pipeline ######


####  Set WD ####
setwd() #Use filepath

#### Load Data and Packages ####
library(tidyverse)
library(sandwich)
library(clubSandwich)
library(dplyr)
library(ggplot2)
library(dotwhisker)
library(scales)
library(devtools)
library(ggeffects)
library(estimatr)
library(fixest)
library(tidylog)
library(patchwork)

load("CredibilityJOP_CleanData.Rda")


####  Additional Data Pre-processing ####


#Lump fatigue, stress, and illness into single variable as pre-registered (presented to respondents separately but theoretically similar)
dat$iv_FSI <- dat$iv_fatigue + dat$iv_stress + dat$iv_illness
dat <- dat %>%
  mutate(
    iv_FSI = if_else(!is.na(iv_FSI), rescale(iv_FSI, to = c(0, 1), from = range(iv_FSI, na.rm = TRUE)), NA)
  )

#### Analyses ####


### Estimate Component Effects ###
# OLS with respondent fixed effects and cluster robust SE

library(fixest)
fit.1 <- feols(dv_eval ~ iv_age + iv_academics + iv_anger + iv_foreignexp + iv_FSI + iv_leaderexp + iv_militaryexp + iv_intel + iv_milstrength | ResponseId, 
               data = dat, vcov = ~ResponseId)

summary(fit.1)


### Estimate Index Additive Effect ###

# Construct indices indicating latent ability and environment 
dat$ability <- dat$iv_age + dat$iv_academics + dat$iv_foreignexp + dat$iv_leaderexp + dat$iv_militaryexp
dat$environment <- dat$iv_intel + dat$iv_FSI

# Rescale the indices to the range 0 to 1, leaving NAs as they are
dat <- dat %>%
  mutate(
    ability = if_else(!is.na(ability), rescale(ability, to = c(0, 1), from = range(ability, na.rm = TRUE)), NA)
  ) %>%
  mutate(
    environment = if_else(!is.na(environment), rescale(environment, to = c(0, 1), from = range(environment, na.rm = TRUE)), NA))


# Estimate additive effect.  OLS with respondent fe and cluster robust SE
fit.2 <- feols(dv_eval ~ ability + environment + iv_anger + iv_milstrength | ResponseId, 
               data = dat, 
               vcov = ~ResponseId)


summary(fit.2)

### Estimate Index Interaction Effect (Test Interactionist Hypothesis)  ###
# OLS with FE, cluster robust SE, and interaction term

fit.3 <- feols(dv_eval ~ ability * environment + iv_anger + iv_milstrength | ResponseId, 
               data = dat, 
               vcov = ~ResponseId)

summary(fit.3)



### Mechanism Tests - Leader Skill  ###
# OLS with respondent fixed effects and cluster robust SE

fit.4 <- feols(dv_skill ~ iv_age + iv_academics + iv_anger + iv_foreignexp + iv_leaderexp + iv_militaryexp | ResponseId, 
               data = dat, 
               vcov = ~ResponseId)

summary(fit.4)

### Mechanism Tests - Situational Difficulty  ###
# OLS with respondent fixed effects and cluster robust SE

fit.5 <- feols(dv_difficulty ~ iv_intel + iv_FSI + iv_milstrength | ResponseId, 
               data = dat, 
               vcov = ~ResponseId)

summary(fit.5)


##### Visual Presentations ######

### Plot Main Component Effects ###

two_brackets <- list(c("Environmental", "Intelligence", "Fatigue / Stress / Illness"),
                     c("Dispositional", "Age", "Academics"))

plot1 <- dwplot(fit.1,
                 vline = geom_vline(
                   xintercept = 0,
                   color = "grey5",
                   linetype = 2
                 ),
                 vars_order = c("iv_age", "iv_leaderexp", "iv_foreignexp", "iv_militaryexp", "iv_academics", "iv_intel", "iv_FSI"),
) %>%
  relabel_predictors(c(
    iv_age = "Age", 
    iv_leaderexp = "Time in Office", 
    iv_foreignexp = "Foreign Policy Experience", 
    iv_militaryexp = "Military Experience", 
    iv_academics = "Academics",
    iv_intel = "Intelligence",
    iv_FSI = "Fatigue / Stress / Illness"
  )
  ) + 
  ggtitle("") +
  xlab("Coefficient Estimate") +
  scale_color_grey() +
  theme_bw()

{plot1 + theme(plot.title = element_text(face = "bold", 
                                         hjust = 1, 
                                         vjust = 3,
                                         size = 14),
               axis.title = element_text(size = 14),
               axis.text = element_text(size = 14),
               legend.position = "none",
               )} %>% 
  add_brackets(two_brackets, fontSize = 1.1)

#Save figure
ggsave("Main Effects Figure.png", width = 8, height = 6, dpi = 300)



### Plot Index Effects ###

## Coefficient plot ##

plot2 <- dwplot(list(fit.2, fit.3),
                vline = geom_vline(
                  xintercept = 0,
                  color = "grey5",
                  linetype = 2
                ),
                vars_order = c("ability", "environment", "ability:environment"),
) %>%
  relabel_predictors(c(
    ability = "Leader Ability", 
    environment = "Decision Environment",
    "ability:environment" = "Ability X Environment"
  )
  ) + 
  ggtitle("") +
  xlab("Coefficient Estimate") +
  scale_color_grey(start = .6,
                   end = .1,
                   labels = c("Additive", "Interactive")) +
  theme_bw()

plot2 + theme(plot.title = element_text(face = "bold", 
                                        hjust = 0.5, 
                                        vjust = 2, 
                                        size = 14), 
              legend.title = element_blank(),
              axis.title = element_text(size = 14),
              axis.text = element_text(size = 14),
              legend.text = element_text(size = 14)
              )

#Save figure
ggsave("Interaction Effects Figure.png", width = 8, height = 6, dpi = 300)


#### Plot Skill and Difficulty Evaluation Models ####

plot3 <- dwplot(list(fit.4),
                vline = geom_vline(
                  xintercept = 0,
                  color = "grey5",
                  linetype = 2
                ),
                vars_order = c("iv_age", "iv_leaderexp", "iv_foreignexp", "iv_militaryexp", "iv_academics", "iv_anger"),
) %>%
  relabel_predictors(c(
    iv_age = "Age", 
    iv_leaderexp = "Time in Office", 
    iv_foreignexp = "Foreign Policy Experience", 
    iv_militaryexp = "Military Experience", 
    iv_academics = "Academics",
    iv_anger = "Anger"
  )) + 
  ggtitle("") +
  xlab("Coefficient Estimate on Leader Skill") +
  scale_color_grey() +
  theme_bw()

plot3 <- plot3 + theme(plot.title = element_text(face = "bold", 
                                        hjust = 0.5, 
                                        vjust = 2, 
                                        size = 14), 
              legend.title = element_blank(),
              legend.position = "none",
              axis.title = element_text(size = 14),
              axis.text = element_text(size = 14),
              legend.text = element_text(size = 14)
)


plot4 <- dwplot(list(fit.5),
                vline = geom_vline(
                  xintercept = 0,
                  color = "grey5",
                  linetype = 2
                ),
                vars_order = c("iv_intel", "iv_FSI", "iv_milstrength"),
) %>%
  relabel_predictors(c(
    iv_intel = "Intelligence",
    iv_FSI = "Fatigue / Stress / Illness",
    iv_milstrength = "Military Strength"
  )) + 
  ggtitle("") +
  xlab("Coefficient Estimate on Environmental Difficulty") +
  scale_color_grey() +
  theme_bw()

plot4 <- plot4 + theme(plot.title = element_text(face = "bold", 
                                                 hjust = 0.5, 
                                                 vjust = 2, 
                                                 size = 14), 
                       legend.title = element_blank(),
                       legend.position = "none",
                       axis.title = element_text(size = 14),
                       axis.text = element_text(size = 14),
                       legend.text = element_text(size = 14)
)


eval_check <- plot3 + plot4 + plot_layout(ncol = 1)

eval_check

ggsave("Skill and Difficulty Evaluation Figure.png", width = 8, height = 6, dpi = 300)


##### Report Results in Tables in Appendix #####

library(modelsummary)
library(flextable)
library(webshot)
webshot::install_phantomjs()

# Dispositional and Environmental Factors

model1 <- list()
model1[["Credibility"]] <- fit.1


tab1 <- modelsummary::modelsummary(model1,
                           coef_map = c(    "iv_age" = "Age", 
                                               "iv_leaderexp" = "Time in Office", 
                                               "iv_foreignexp" = "Foreign Policy Experience", 
                                               "iv_militaryexp" = "Military Experience", 
                                               "iv_academics" = "Academics",
                                               "iv_intel" = "Intelligence",
                                               "iv_FSI" = "Fatigue / Stress / Illness",
                                              "iv_anger" = "Anger",
                                               "iv_milstrength" = "Relative Military Strength"),
                           gof_omit = 'AIC|BIC|Std.Errors|R2 Adj.',
                           notes = "OLS estimated with respondent-fixed effects and cluster-robust SE.",
                           output = 'flextable'
                           )

tab1 <- tab1 %>% 
  autofit() %>%
  #add_header_row(values = c("Effects of Dispositional and Environmental Factors on Credibility"), top = TRUE, colwidths = 2) %>%
  font(fontname = "Helvetica")


# Indices

models2_3 <- list()
models2_3[["Additive"]] <- fit.2
models2_3[["Interactive"]] <- fit.3

tab2 <- modelsummary::modelsummary(models2_3,
                                   coef_map = c(    "ability" = "Leader Ability", 
                                                    "environment" = "Decision Environment", 
                                                    "ability:environment" = "Ability x Environment"),
                                   gof_omit = 'AIC|BIC|Std.Errors|R2 Adj.',
                                   notes = "OLS estimated with respondent-fixed effects and cluster-robust SE.",
                                   output = 'flextable'
)

tab2 <- tab2 %>% 
  autofit() %>%
  #add_header_row(values = c("Effects of Dispositional and Environmental Factors on Credibility"), top = TRUE, colwidths = 2) %>%
  font(fontname = "Helvetica")

# Mechanism Probe

model4_5 <- list()
model4_5[["Skill"]] <- fit.4
model4_5[["Difficulty"]] <- fit.5

tab3 <- modelsummary::modelsummary(model4_5,
                                   coef_map = c(    "iv_age" = "Age", 
                                                    "iv_leaderexp" = "Time in Office", 
                                                    "iv_foreignexp" = "Foreign Policy Experience", 
                                                    "iv_militaryexp" = "Military Experience", 
                                                    "iv_academics" = "Academics",
                                                    "iv_anger" = "Anger",
                                                    "iv_intel" = "Intelligence",
                                                    "iv_FSI" = "Fatigue / Stress / Illness",
                                                    "iv_milstrength" = "Relative Military Strength"),
                                   gof_omit = 'AIC|BIC|Std.Errors|R2 Adj.',
                                   notes = "OLS estimated with respondent-fixed effects and cluster-robust SE.",
                                   output = 'flextable'
)

tab3 <- tab3 %>% 
  autofit() %>%
  #add_header_row(values = c("Effects of Dispositional and Environmental Factors on Credibility"), top = TRUE, colwidths = 2) %>%
  font(fontname = "Helvetica")

# Save tables

save_as_image(tab1, path = "~tableA.1.png")  #Insert full working directory
save_as_image(tab2, path = "~tableA.2.png")
save_as_image(tab3, path = "~tableA.3.png")
